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We extend earlier work [Phys. Rev. Lett. 84, 3740 (2000)] on the statistical mechanics of the 
cubic one-dimensional discrete nonlinear Schrodinger (DNLS) equation to a more general class of 
models, including higher dimensionalities and nonlinearities of arbitrary degree. These extensions 
are physically motivated by the desire to describe situations with an excitation threshold for creation 
of localized excitations, as well as by recent work suggesting non-cubic DNLS models to describe 
Bose-Einstein condensates in deep optical lattices, taking into account the effective condensate 
dimensionality. Considering ensembles of initial conditions with given values of the two conserved 
quantities, norm and Hamiltonian, we calculate analytically the boundary of the 'normal' Gibbsian 
regime corresponding to infinite temperature, and perform numerical simulations to illuminate the 
nature of the localization dynamics outside this regime for various cases. Furthermore, we show 
quantitatively how this DNLS localization transition manifests itself for small-amplitude oscillations 
in generic Klein-Gordon lattices of weakly coupled anharmonic oscillators (in which energy is the 
only conserved quantity), and determine conditions for existence of persistent energy localization 
over large time scales. 

PACS numbers: 05.45.-a, 63.20.Pw, 63.70.+h, 03.75.Lm 



I. INTRODUCTION 

There is a large interest in many branches of current science in the topic of localization and energy transfer in 
Hamiltonian nonlinear lattice systems (see e.g. Ref. for a comprehensive review, and Refs. 0- 13- DJ for more recent 
progress). Under quite general conditions, such lattices sustain exact, spatially exponentially localized and time- 
periodic, solutions termed intrinsically localized modes (ILMs) or discrete breathers (DBs). Although their existence 
as exact solutions has been rigorously proven in many explicit cases (|3, |g, and e.g. Ref. pi |a |9| and references 
therein for extensions), there is still an ongoing debate regarding their relevance to actual physical phenomena, at 
nonzero temperatures. Important fundamental questions concern whether ILMs may exist in thermal equilibrium, 
or if not, whether their typical lifetimes are long enough to considerably influence transport properties of crystals, 
biomolecules, etc. 

A frequently studied example of a non-integrable Hamiltonian lattice model is the discrete nonlinear Schrodinger 
(DNLS) equation (see Refs. [nj El for recent reviews of its history, properties and applications). This model is of 
great interest from a general nonlinear dynamics point of view, where it provides a particularly simple system to 
analyze fundamental phenomena such as energy localization, wave instabilities etc., resulting from competition of 
nonlinearity and discreteness, as well as from a more applied viewpoint describing e.g. arrays of nonlinear optical 
waveguides or Bose-Einstein condensates in external periodic potentials. The DNLS equation can be derived through 
an expansion on multiple time-scales of small-amplitude oscillations in a generic class of weakly coupled anharmonic 
oscillators [Klein-Gordon (KG) lattice], and thus approximates the KG dynamics over large but finite time-ranges 
(see, e.g., Refs. 0,0,^3)- A particular feature of the DNLS model is the existence of a second conserved quantity 
in addition to the Hamiltonian: the total excitation number (norm) of the solution. In the KG model, this quantity 
roughly corresponds to the total action integral, which thus must be an approximate invariant in cases where the 
DNLS description of the KG dynamics is acceptable. 

A fundamental question is, for which kinds of spatially extended initial states may we expect spontaneous formation 
of persistent localized modes in such lattices? The answer generally requires a statistical-mechanics description of the 
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model. Due to the existence of a second conserved quantity, it has been possible to obtain some analytical results for 
the thermodynamic properties of the DNLS model in the grand-canonical ensemble, by identifying the norm with the 
number of particles in the standard Gibbsian approach (this is also its physically relevant interpretation in the Bose- 
Einstein DNLS realization). In Ref. 1.4], it was found that the onset of persistent localization could be identified with 
a phase transition line in parameter space, such that on one side the system thermalized according to the Gibbsian 
distribution with well-defined chemical potential and (positive) temperature, while on the other side the dynamics 
was associated with a negative-temperature behavior (for finite systems) creating a small number of large-amplitude, 
standing localized breathers. The transition line was shown to correspond to the limit of infinite temperatures in the 
'normal' regime. Similar properties were later found also for other types of lattice models with two conserved quantities 
in Ref. 01 • Most recently, in Ref. ^(| Rumpf revisited the statistical-mechanics description of the DNLS localization 
transition. Under the particular assumption of small-amplitude initial conditions (an assumption not made in Ref. 
|l4j). he argues that the phase space generally can be divided into two weakly interacting domains, corresponding to 
low-amplitude fluctuations ('phonons') and high-amplitude peaks (breathers), respectively. Explicit expressions for 
macroscopic quantities, valid not only in the 'normal' regime but in the full range of parameter space, can then be 
obtained by assuming the two domains to be in thermal equilibrium with each other, and the emergence of localized 
peaks in the 'anomalous' phase arises as the system strives for maximizing its total entropy. Under these conditions, 
the temperature is not negative but infinite in the thermal equilibrium state with coexisting large-amplitude breathers 
and small- amplitude fluctuations. 

Let us mention a number of reasons that have lead us to revisit and extend the results of Ref. 14:]. First, so far only 
the one-dimensional (ID) case with cubic nonlinearity was considered. However, apart from the natural interest in 
considering two- and three-dimensional (2D and 3D) physical situations, there is also a fundamental difference to the 
ID case: there is an excitation threshold for creation of localized excitations for the cubic DNLS model [l7i fill IT^ | . A 
similar threshold also occurs in the ID DNLS equation for noncubic nonlinearities of the form \ip m \ 2a with a > 2 
[TtI fl8L Il9l | , and generally the condition aD > 2 for existence of an excitation threshold in D dimensions is the same 
as the condition for collapse of the ground-state solution of the corresponding continuous NLS equation (e.g. Ref. 
ppf). For this reason, one sometimes studies the ID DNLS equation with larger a hoping to capture the main effects 
of higher dimensionality in a simpler ID model (e.g. Refs. |2ll l22l l23l |24|L Recently |2jj, similar arguments were 
also used in the study of a ID KG chain with a <fi 8 on-site potential, to mimic the effects of an excitation threshold 
for breathers in the thermalization dynamics of a three-dimensional KG-lattice. (A similar relation between degree 
of nonlinearity and dimension is valid also for KG-lattices, see Ref. an d R^f- US f° r recent extensions.) Thus, 
it is of interest to investigate the nature of the statistical localization transition for various degrees of nonlinearity 
and dimensions, in order to elucidate (i) whether it is qualitatively affected by the existence of a breather excitation 
threshold, and (ii) whether quantitative effects arising from increasing a agree with those from increasing D. 

While the above connection motivates the study of particular on-site nonlinearities with a = 2 and a = 3, recent 
progress in studies of Bose-Einstein condensates in optical lattices also provide motivation for considering non-integer 
values of a < 1. It has namely been suggested [27|], that the effective power of the nonlinearity in the tight-binding 
DNLS approximation depends on the effective dimensionality d of the condensate in each well, such that a — 2/(2 + d) 
where d = Q, 1, 2, or 3. Moreover, it is tempting to suggest a connection between the statistical localization transition 
in the DNLS model and experimentally observed superfluid-insulator transitions of the condensate (e.g. Refs. psll^ 
and references therein). 

Last, but not least, we wish to employ the results for the DNLS model to give quantitative predictions for breather 
formation in generic KG models, and in particular describe what kinds of initial conditions yield long-lived breathers in 
the regime of weak coupling and small averaged energy density where the DNLS approximation is justified. Although 
particular examples of the manifestation of the DNLS localization transition in KG models have been given earlier 
|13| , we here derive explicit general approximate expressions for the transition line in terms of direct properties of the 
KG initial state. Due to the violation of norm (or action) conservation, the transition in the KG model is not strict, 
and we perform numerical simulations to investigate how the long-time dynamics is influenced by the slow variation 
of the almost conserved quantity. We suggest that the approach proposed here could be used to clarify the findings 
regarding the role of breathers in thermalized KG lattices (with or without energy gaps) of Refs. |25l l30| . which did 
not employ the connection to the DNLS model. 

The structure of this paper is as follows. Sec. [H] describes the statistical mechanics of general DNLS models. Sec. 
Ill Al generalizes the statistical- mechanics approach of Ref. 14] to ID models with general degrees of nonlinearity. As 
particular examples, we consider initial conditions taken as traveling (Sec. Ill A 1|) and standing (Sec. Ill A 2|) waves. 
We obtain simple analytical conditions for the transition into the statistical localization regime, and illustrate with 
numerical simulations the actual dynamics on both sides of the transition. Sec. Ill Bl extends these results to higher 
dimensions. In Sec. IIIII we describe how the results from the DNLS model can be transfered into approximate 
conditions for statistical formation of long-lived breathers in weakly coupled Klein-Gordon chains, and confirm and 
illuminate these predictions with numerical simulations. Sec. lIVI gives some concluding remarks and perspectives. 
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II. STATISTICAL MECHANICS OF GENERAL DNLS MODELS 



A. ID model with general degree of nonlinearity 

Generalizing the ID DNLS equation of Ref. to include a nonlinearity of arbitrary (homogeneous) degree, we 
consider the DNLS equation in the form: 

ilp m + C(tp m+1 + 1p m -l) + \lpm\ 2rT 1pm = 0, (1) 



and norm (exci- 



with the two conserved quantities Hamiltonian H — J2 m CCV'mV'm+i + V'mVWfi) + ~^fi\i J m\ 2 ' J+2 

tation number) A = Ylm I^Aml 2 - Compared to Eq. (1) of Ref. [T3| . we have used v = 1 as the coefficient of the nonlinear 
term, included a coupling constant C > in front of the coupling terms, and generalized \ip m \ 2 i)m to \ip m \ 2(T ipm with 
cr > 0. Note that although we formally discuss the case of positive intersite coupling and positive nonlinearity, this is 
not a restriction, since changing the sign of C is equivalent to the transformation ip m — > (— l) m ip m , while the same 
transformation followed by a time reversal t — s- — t is equivalent to changing the sign of the nonlinearity. Thus, all 
obtained results can be directly transformed to the cases of negative coupling and/or negative nonlinearity. Any finite 
coefficient in front of the nonlinear term can also be obtained through a simple rescaling. 

With a canonical transformation into action-angle variables, ip m — y 1 A m e %l ^ m , the Hamiltonian for a chain of N 
sites becomes 

N / i \ 

H = (2CVA m A m+1 cos(0 m - cf* m+1 ) + — -A^ 1 J , (2) 

m— 1 ^ ' 

and the norm 

AT 

A=Y^ A m- (3) 



m— 1 



We first note that the staggered (q = it) stationary homogeneous plane-wave solution V£r™ } = ^/A/Ne ml *e iAt , 
with A = -2C + {A/NY, minimizes H at fixed A and N, for all a. The minimum value is thus H {min) = -2CA + 
-^A a+1 /N a . To prove this, write 



N 



n - n {min) = Y 



2C (^/A m A m+1 cos(0 m - <j> m+1 ) + + — j-j- I A^+ l - (jj 



The first part is positive, since {\J A rn A m+1 cos(</> m - (f> m +i) + 4) ^ S (4 ~ y/A m A m+1 ) = 
i ^2 (V A m — \/A m+ i) > 0. The second part is also positive, which can be seen from Holder's inequality: 

EMfcl < (E\ak\ p ) 1/P (E\h\ q ) 1/q , if 1/p+l/q = 1. Let a fc = A m A = l,p = tr + l.g = 1 + 1/(7, which gives 
^A^ 1 " 1 — (E J 4 m ) <T+1 > 0. Notice also that "ft( mm ) is bounded from below as a function of A for any finite 
number of sites N, with the global minimum H {min) = --^- [ N{2C) 1+1 / a obtained for A = N{2C) 1 / a . 

Similarly to the work of Ref. |l4j , we use standard Gibbsian statistical mechanics to predict macroscopic average 
values in the thermodynamic limit, by treating the norm A as analogous to 'number of particles' in the grand-canonical 
ensemble. As in Eq. (2) in Ref. [l4( , the grand-canonical partition function is thus defined as 

OO /-27T N 

[] dcb m dA m e-f^ A \ (4) 

771—1 

where f3 = 1/T (in units of ks = 1) and play the roles of inverse temperature and chemical potential, respectively. 
Using (0-® an d integrating over the phase variables (j) m yields 

roc / a& \ 

Z = (2ir) N / T\dA m I {2[3C^ A m A rn+1 )e- pAm y^ t+ ' 1 ) , (5) 
J o m 

where Io(z) = i f£ e zcos0 d6 is the modified Bessel function of the first kind. From this expression, one could proceed 
as in Ref. |l4j by symmetrizing the partition function and using the transfer integral operator to obtain thermodynamic 
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quantities in the limit N — > oo, corresponding to the regime in (A,7i) parameter space with well-defined chemical 
potential and (positive) temperature. This is however not our main purpose here. Instead, we focus on the phase 
transition line defined by the boundary of this regime (/? = 0, u_= oo, with /3/j, = 7 finite), which signals the transition 
into the regime of persistent localization, suggested in Ref. |l J | to be associated with a negative-temperature type 
behavior for finite lattices and time-scales. 

Close to the high-temperature limit (3 — ► + , we can approximate the slowly increasing Bessel function with Iq sw 1 
(which is mathematically equivalent to letting C — > 0, corresponding physically to thermalized independent units). 
The partition function then becomes Z ~ (2iry(f3, [i)) N , where 



(3x a+l 1 ( (3x a+1 



fa v + lj dX+ 2(a + l)*J 



(7+1 2 \ (7+1 

2 



dx 



x 



But / °° x n e ax dx — r ^"^ 1 1 - > (where T is the Gamma function, T(n + 1) = n! for integer n). This yields 

1_ (3 I> + 2) 1 /3 2 r(2a + 3) 
V[P ' N fa a + l {fa) a + 2 2(a + l) 2 (/3 M ) 2ff+3 "' 

Thus, close to the limit of (3 ~ * 0, /1 — ► 00 with /3/x = 7 constant, we can neglect all higher-order terms in /3, and 

JY 



obtain y{j3,fj) ~ 4 fa ^ w+2 ■ Finally, for the partition function in the high-temperature limit we get 



1 / pr{a + i) . 
(2 " } (^fI^W^' ■ (6) 



For small /3 this reduces to InZ ~ 7Vln(27r) — A r ln(/?//) — N ^^tll , so that we have in the high-temperature limit 
for the average energy: 



and for the average norm 



c A > = IdlnZ^N NT(a + 2) 

13 dfi ~ Pn ^{fay +1 ' { ' 

(The second term here is negligible.) Thus, the relation between the energy density h = < ^ > and the norm density 
a = < ^ > in the high-temperature limit is 

h = r((7 + l)a a+1 , (9) 

where T(cr + 1) can be replaced by a\ for integer a. Note that the quantity 7 = (3fi indeed is well-defined and finite 
in the high-temperature limit for any nonzero norm density, 7 ~ 1/a according to ||SJ). 

For a = 1, the corresponding phase diagram was illustrated in Fig. 1 of Ref . [lj. Thus, for any given norm density 
a, typical initial conditions with (Hamiltonian) energy density h smaller than the critical value © are expected to 
thcrmalize (after 'sufficiently' long times) according to a Gibbsian equilibrium distribution at temperature T = 1/(3 
and chemical potential \x. The correspondence between (a, h) and (f3, //) generally has to be found numerically through 
the transfer integral formalism as in Ref. [T^j . but in the small-amplitude limit a — * analytic expressions can be 
obtained as shown in Ref. [16|, Eqs. (7)- (8). Numerical evidence that such a thermalization generally takes place after 
sufficiently long integration times was given in Fig. 2 of Ref. 14] (for (7 = 1). 

On the other hand, for initial conditions with energy density h larger than the critical value © this description 
breaks down, and one finds numerically that persistent large-amplitude standing breathers are created. Heuristically, 
this can be understood as follows: For fixed norm A, it is generally possible to maximize the Hamiltonian H, and 
the maximizing solution is a single-site peaked, exponentially localized stationary standing breather (see, e.g., Ref. 
fifty), which for large A becomes essentially localized at one site so that J{( max ) ~ A a+1 /(cr +1). Considering in the 
microcanonical ensemble (fixed A, Tt and N) the entropy S(Tt , A, N) (i.e. the logarithm of the number of microstates) 
as a function of Tt, it is zero at Tt^" 1171 ^ (A, N) defined above, increases towards its maximum when © is fulfilled and 
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T = oo (since 1/T = 8S/&H\a,n), and then again decreases towards zero at TC^ max \A). Thus, in the microcanonical 
ensemble at finite A and N, the temperature is well-defined and becomes negative when h = Ti/N is larger than the 
critical value 10. Returning to the grand-canonical ensemble, it is then possible for the part of the system which is in 
the negative-temperature regime to increase its entropy by transferring some of its superfluous energy into localized 
breathers, which consume only a small amount of the norm. In other words, the 'overheated' negative-temperature 
system 'cools itself off' by creating breathers as 'hot spots' of localized energy. Such a mechanism for energy localization 
works quite generally in systems with two conserved quantities (see, e.g., Ref. |15| and references therein). Indeed, 
this type of argument could be used to explicitly calculate the thermodynamic properties of the DNLS model in the 
limit of small a. where phase space naturally divides into a small-amplitude 'fluctuation' part and a large-amplitude 
'breather' part |l6l | which only interact weakly. In that case, the equilibrium state which maximizes the total entropy 
for h larger than the critical value 10 should consist of one single breather, with the rest of the lattice corresponding 
to an ordinary Gibbsian distribution at T — oo fl6l | (although the numerical simulations in Ref. never reached 
such a state, but rather one with a finite breather density). However, when a increases, the large-amplitude and 
small-amplitude parts will not separate straightforwardly anymore, and the thermodynamic equilibrium properties 
for general a remain unknown. Some of the numerical simulations reported below aim at shedding some light on this 
issue. 

Let us now discuss the thermodynamical equilibrium distributions for some particularly interesting choices of initial 
conditions. For certain families of exact solutions we can analytically compute the curves h (a), and thus within these 
families obtain the transition into the statistical localization regime by finding their intersections with the phase- 
transition line I0- Evidently, if the initial condition is strictly an exact solution thermalization will not occur, but 
often solutions are linearly unstable e.g. through modulational ^3 or oscillatory [l3Ll3lLl3^| instabilities, which may 
cause rather rapid thermalization (see e.g. examples for a = 1 in Refs. [LH. l3l|). Even for weakly perturbed linearly 
stable solutions as initial conditions, it is expected that generically nonlinear instability mechanisms finally should 
lead to thermodynamic equilibrium; however the equilibration times can be extremely long as Arnol'd-type diffusion 
processes are involved. 

In the numerical investigations below, we mainly focus on the distribution function p(A m ) for the amplitudes 
A m — \ijj m \ 2 , which most clearly illustrates the localization properties. In the standard Gibbsian regime, the statistical 
prediction for p(A m ) can also be obtained through the transfer integral formalism as show in Ref. [l4j- Here, let us 
only note that close to the high-temperature limit (3 — > 0, this prediction yields (again by approximating Iq ps 1) 

logp(A m ) ~ -jA m - f3AH 1 /{a + 1), (10) 

i.e., the curvature is zero for (3 = and becomes negative (positive) for positive (negative) temperatures. Thus, 
negative temperatures favor large-amplitude excitations. 



1. Traveling waves 

For a traveling wave, which is an exact solution of the form ip m — ^/ae iqm e (with A = 2C cose/ + a a ), we have 

h = 2Cacosq H . (11) 

(7+1 

Similarly as for the well known case a = 1 [lj, traveling waves with |g| < tt]2 are modulationally unstable and those 
with 7r/2 < \q\ < tt linearly stable also for general a > (see, e.g., Ref. |22|). To find when such a solution crosses 
the (3 = curve we put l(TT)l equal to ©, which yields 

2(a + l)Ccos g 

r(a + 2)-i ' v 1 

where, as before, r(<r + 2) can be replaced by (cr + 1)! for integer a. Thus, for any a and \q\ < tt/2, there is a threshold 
value for the norm density given by (|12|l . so that only above this threshold, one will be in the 'normal' Gibbsian 
positive-temperature regime, while below it we expect statistical localization. The predicted threshold, plotted in Fig. 
^ becomes quite small for large cr due to the factorial in the denominator, but increases rapidly for a smaller than 1 
(e.g. for cr = 0.4, corresponding to 3-dimensional Bosc-Einstein condensates in the model of Ref. [27]], the threshold 
is a ps 455 for q — and C — 1). On the other hand, for tt/2 < \q\ < tt, one is always in the normal thermalizing 
regime. 

In Fig. |3 we show some examples of resulting distribution functions p{A) obtained from long-time numerical inte- 
grations of constant-amplitude (q = 0) initial conditions. For the small value a = 0.4 (Fig. 2(a)), we can note that the 



6 



500 



400 



300 



200 



100 




0.2 



0.4 



0.6 



0.8 

q 



1.2 



1.4 



FIG. 1: Maximum value of norm density a for statistical localization to occur, according to 1121 1. from initial condition being 
a travelling wave of wave vector q, for ID DNLS with various a. From top to bottom at q — 0: a = 2/5, 1/2, 2/3, 1, 2, 3. Inset 
is blow-up for small a. 




FIG. 2: Numerically obtained distribution functions p(A m ) resulting from long-time integration (t — 1.1 • 10 6 ) of (unstable) 
constant-amplitude initial states ipm,(0) = •/a, for a ID DNLS chain with N — 10000 (C = 1). (a) a — 0.4, a — 100 (squares), 
a — 455 (triangles), and a — 1000 (circles), (b) a — 3, a — 0.6 (circles), and a = 0.8 (squares). The dashed lines represent best 
fits to ftfy. 



numerics perfectly confirms the predicted transition at a ~ 455 (with a linear dependence log p(A) ~ — according 
to Hl()(l ). However, to achieve an appreciable difference between the distributions at either side of the transition point 
(compared e.g. to the case a = 1 illustrated in Figs. 2-3 in Ref. |l4|]), we had to choose initial conditions quite far from 
the transition line. Then, fitting 7 and /3 in Eq. i|10|) to the obtained distributions, we find small values of j3 with the 
expected (opposite) signs in the two cases. We attribute the smallness of j3 even for values of a far from the transition 
point to the weakness of the nonlinear effects for small a. Moreover, as we illustrate with another example in the 
following subsection, the thermalizing dynamics in the localization regime is extremely slow for small a. Although 
we can clearly identify several breather-like excitations with amplitudes considerably higher than their surroundings 
in the simulations for a < 455, they are generally not persistent but transient and recurring. Thus, it is necessary 
to remember, that curves such as those for a < 455 in Fig. GJa), obtained after long but finite-time integrations, do 
generally not represent true equilibrium distributions in the localization regime, but rather an intermediate stage in 
the approach to equilibrium by breather-forming processes in a negative-temperature regime. The a = 3 case (see Fig. 
|2Ib)) contrasts this by showing an appreciable number of persistent breather- like excitations in the breather forming 
regime a = 0.6 (circles). For a = 3 the critical amplitude is a ~ 0.7, and we see that the distribution functions obey 
the predicted behavior Eq. I|10|l both in the breather forming and in the normal regime (a = 0.8) (squares) until finite 
size effects set in at A ~ 4. 
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2. Standing waves 

In addition to travelling waves, there are also exact solutions in the form of standing waves (SWs), which are 
time-periodic non-propagating (i.e., with their complex phase spatially constant) solutions, with an inhomogeneous 
amplitude distribution \ip m \ 2 being periodic or quasiperiodic in space [ill. l3ll l32| . In the linear limit a — > 0, a standing 
wave of wave vector Q (0 < \Q\ < n) is a linear combination of two counterpropagating travelling waves q = ±Q, 
i.e., tp m ~ \f2asm-(Qm + ip)e lAt for small a. As a increases, one finds 0,0, H2, that only for particular phases tp 
can these linear SWs be continued into exact nonlinear SW solutions. These can be divided into two distinct classes: 
phases <p = ±(ir — Q)/2 — m'Q (to' integer) continue into solutions called 'type E, while either <p — —m'Q (for generic 
Q) or <p = —(rn' + i)Q (for special Q = J^jn n, k,k' integers) yield solutions called 'type IT. (These two types 
of solutions can be represented as elliptic and hyperbolic cycles, respectively, of the cubic real 2D map 0, 0, 
when cr = l). In physical space, they are distinguished by their positioning in the lattice, with type-E SWs centered 
symmetrically between lattice sites at to = m'+ 1, and type-H SWs centered antisymmetrically either around a lattice 
site at to = to' (generic Q) or between sites at to = ml + | (for Q = ^txjTt), respectively. In the opposite limit of 
large a, which is mathematically equivalent to C — > 0, both classes of solutions can be generated from a circle map, 
distributing solutions if> m = Q,±\fAe lA " t periodically or quasiperiodically in space \13L l3ll l32| . Type-E solutions are 
generally linearly unstable, while type-H solutions are linearly stable for large a/C but generally oscillatorily unstable 
for small a/C (for a = 1, see Refs. [13, HH E3 ) • 

Particularly interesting in this context are the SWs with Q = tt/2, which have the form ip2n+i = 0, 7/^+2 = 
(-l)"v / 2ae^ 2a ^ t (type H), and -02n+i = V'2n+2 = (-1)" \fae la " * (type E), respectively. For small a, any wave 
(travelling or standing) with wave vector 7r/2 coincides with the phase transition line 10 as noted in Ref. [16j. This 
is not true in general, and in particular it is clear from l|12(l that a travelling wave with q = ir/2 lies inside the regime 
of 'normal' thermalization for all nonzero a. On the other hand, it was noticed in Refs. |ll l31| . that for a = 1 the 
curve h(a) for the Q = tt/2 type-H SW indeed coincides with the phase transition line and that type-H SWs 
with \Q\ < 7r/2 generally resulted in creation of large-amplitude breathers, and those with ir /2 < \Q\ < ir in 'normal' 
thermalization (see e.g. Figs. 6-9 in Ref. |13|). However, for general a we now have the relation for Q = tt/2 type-H 
SWs: 

h = ^-a°+\ (13) 
o + 1 

Thus, only for the particular case a — 1 considered in Refs. [l3l l3l| do the coefficients in and (|13|) agree. In 
general are the transition line into the phase of statistical localization and the line defined by the Q = ir/2 type-H 
SW different. For < a < 1 the ir/2 type-H standing wave will always be in the breather-forming regime, while for 
a > 1 it will always be in the normal thermalizing regime. This is illustrated by the numerical simulations in Fig. EJ 

For cr = 3, Fig. [3^a) clearly confirms a positive-temperature behavior, with a distribution function well fitted 
by HlOf) with positive (3, and very small probability for large- amplitude excitations. For a = 2/3 we do observe, as 
predicted, a small positive curvature of the distribution function at finite times, as well as a tendency towards creation 
of large-amplitude breathers (e.g. the four points between A = 12 and A = 14 in Fig. EJa)). However, even for very 
large systems and long integration times, the breathers found are not persistent but transient and recurring, as for 
the small-cr case discussed in the previous subsection. 

To check to what extent the finite-time averaged distribution functions in Fig. Ola) are representative for the true 
equilibrium distributions, we monitor the average of the contribution to the total Hamiltonian from the coupling 
part, h coup (first term in Eq. ©). By definition, < h coup >= in equilibrium at the transition line (3 — 0, and, by 
the particular choice of n/2 type-H SWs as initial conditions, h coup (0) = for all a. For (7 = 1, Fig. |3{b) confirms 
that < /icoup >, although being positive for intermediate times, asymptotically approaches zero as expected. For 
ct = 3, < /i C oup > approaches asymptotically a negative value, which is typical in the positive-temperature regime, 
and implies a preference for out-of phase excitations at neighboring sites. For a = 2/3, a superficial look at the 
main Fig. EJb) seems to indicate an asymptotic approach to a strictly positive < h coup >, signifying a preference 
for in-phase excitations at neighboring sites. However, as is shown by the inset in Fig. |3fb), the simulation indeed 
has not reached a stationary regime even after t = 2 • 10 5 , and there is a very slow decrease, close to logarithmic 
in time, of < h coup >. We attribute this to an on-going process of formation of large-amplitude breathers. Note 
that, if the hypothesis of approaching a thermodynamic equilibrium state consisting of one (or a finite number of) 
breather(s) together with an infinite-temperature phonon bath would be correct, we should always asymptotically 
have < h coup >= in the breather- forming regime for N — > oo. Thus, our simulations are consistent with (although 
by no means proving) this hypothesis. However, extrapolating the tendency of the curve in the inset in Fig. ffib) 
to larger times would yield < h coup >= only after t ~ 10 70 , i.e., the times to reach a true equilibrium state in 
the breather-forming regime are indeed extremely long! Let us only for completeness stress, that the observed slow 
decrease of < h coup > is a true behavior of the system, and not an artifact of numerical drifting of the conserved 
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FIG. 3: (a) Time-averaged (non-normalized) distribution functions p(A m ) for weakly perturbed ty/2 type-H SWs with a = 1 
(..., — v2) 0, V% 0, ...) as initial conditions (C = 1). (*): a — 2/3; (+): a = 1; (x): cr = 3. The points are obtained by averaging 
over 96 time instants for 10000 < t < 200000 and N = 10000 (cr = 2/3), 177 time instants for 500000 < t < 1380000 and 
N = 1000 (<t = 1), and 91 time instants for 10000 < t < 55000 and N = 1000 (a = 3). Straight lines for cr = 2/3 and a = 1 are 
predictions from HUH with /3 = and 7 = 1/a = 1, while curve for a = 3 is prediction from (1101 with fitted values of f3 — 0.042 
and 7 = 0.65. (b) Average (over space and time) of the coupling part of h (i.e., < 2C V A m A m +i cos(c/) m — m +i) >) versus 
time for the simulations in (a) with, from top to bottom at t — 10000, cr = 2/3, a = 1, and cr = 3, respectively. Magnification 
in inset illustrates the slow long-time decrease for a — 2/3. 



quantities during the simulation time. Indeed, there is a slow numerical drift of h (increasing approximately 4 • 10 -12 
per time unit), but this is negligible compared with (and in addition in the opposite direction to) the tendency in Fig. 
OUb) over the used integration time. 

In this context, we should also remark that, in contrast to the ordinary DNLS case a — 1 where the Q = w/2 
type-H SWs are always linearly unstable for small a, this is not the case for < a < 1/2 where they are linearly 
stable for all a. It follows from a standard linear stability analysis (cf. e.g. Ref. H3|), that these solutions (also termed 
'period-doubled states' in Ref. 34|) are oscillatorily unstable for small- wavelength relative perturbations when the 
condition (2a) 2rT + 16(1 — 2a) < is fulfilled, and linearly stable otherwise. Note that this condition is always fulfilled 
for small a if a > 1/2, but can never be fulfilled if cr < 1/2. 

Regarding the type-E SW with Q = tt/2, we note that this solution is a special case of the general class of equivalent 
solutions ■02n+i = (—l) n \Zae la *, i/> 2n.+ 2 = (— V) n \fae taa e %a *, where ao can take any real value (this class of solutions 
were called V — 7r states' in Ref. and 'phase states' in Ref. [13). Putting ao = yields the type-E SW with 
Q = tt/2, while ao = 7r/2 yields the travelling wave with q — tt/2. Thus, h = a a+1 j(a + 1) for all solutions in this 
class, and they belong to the 'normal' thermalizing regime for all nonzero a. 



B. Higher-dimensional models 



An important point to note is, that the results from the previous subsection are readily generalized to higher- 
dimensional DNLS equations. Considering e.g. the 2D case for a quadratic lattice of N sites, we can write the 
expression for the Hamiltonian analogous to (J2J) as 

v^v , 1 . 

'H= \2C y/ A m ^ n A m+1 _ri COs((f) m<n - 4> m +l,n) + \/ A m , n A m , n+ i COs((/) m . n - <j) m ,n+l) + ~ ^ ^mtn f ■ ( 14 ) 

m.n—1 ^ ' 

With this Hamiltonian, the expression for the grand-canonical partition function analogous to (JSJ) becomes 

^00 Vn 



Z = (2ir) N f ]J dA m<n I (2pC^A m>n A m+1<n )I (2pCy/A m!n A m>n+1 )e 0Am ' n {^ + ») } (15 ) 

•* m,n=l 

from which we obtain the behavior close to the high-temperature limit /3 — > + again by approximating Iq w 1. 
Thus, in this limit all results are independent of dimension, which is a consequence of the equivalence of this limit 
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FIG. 4: (Color online) Resulting distribution functions p(A m ) after long-time integrations of initial conditions consisting of 
weakly perturbed 2D constant-amplitude (q x — q y = 0) unstable solutions with (a) a = 5, (b) a — 7, and (c) a = 8. Curves in 
(a)-(c) have been obtained by averaging over a number of different random initial perturbations (8 in (a), 14 in (b), and 100 
in (c)); system size 128 x 128 (N = 16384); integration times t = 500000 (a), 200000 (b), respectively 50000 (c). Curves in (d) 
have been obtained from one single realization for a 50 x 50 system with a = 7, by averaging over 20 different time instants in 
the intervals 500 < t < 10000 (squares), and 110000 < t < 300000 (circles), respectively. The scales are such that, in (a)-(c) 
the dots with smallest probability correspond to one site in one realization, and in (d) to one site at one time instant. Straight 
lines are predictions according to 1101 with (3 = 0. (<r = C = 1.) 



to C — > 0, i.e., thermalized independent units which neglect all interaction terms. Thus the expression @ for the 
phase-transition line is indeed valid for given <r in any dimension! 

To take a specific example in 2D, consider again a travelling plane wave ipm,n = v /ae^ 9x " l+9yn - ) e iAt (with A = 
2C{cosq x + cosq y ) + a a ). It follows from standard analysis (see, e.g. Ref. [23), that the travelling waves are linearly 
stable only if tt/2 < \q x \, \q y \ < tt, and modulationally unstable if either or \q y \ (or both) are smaller than n/2. 
We immediately obtain the expression for the Hamiltonian density by just replacing cos q with cos q x + cos q y in the 
ID expression (JTTJ) , and likewise we obtain the expression for the statistical localization transition analogous to l|12|l : 

a<T = 2(a+l)C(cosq x +cosq y ) 

r(cr + 2)-l 1 ' 

Thus, a necessary condition for breather formation from 2D travelling waves is to have cos^ + cos<?j, > 0, i.e., either 
\Qx\ or \ly\ (but not necessarily both) has to be smaller than 7r/2. Just as for ID, the dynamics always enters the 
'normal' thcrmalizing regime if the norm density is large enough, and the largest possible a for breather formation 
occurs for q x — q y — 0. Note that for this constant-amplitude solution, the threshold in a for a = 1 is multiplied by 
an additional factor of 2 compared to the analogous ID q = case in ljT2jl. becoming 8C instead of 4C. 

Numerical illustrations of the resulting distribution functions in either regimes, together with predictions according 
to JTOJl, are shown in Fig. Note that in the breather-forming regime (Fig. 01(a), (b), and (d)), the distributions 
closely follow the straight lines p{A m ) = ±-e~ Am / a corresponding to /3 = in (|1(J|) up to some threshold value of A m . 
We find, that extending the integration time this breaking point typically moves in the direction of larger A m . For 
small integration times, one finds a smooth curve with positive curvature, indicating a negative-temperature behavior 
as discussed in Ref. [l4|- However, for larger times the tendency is that the curve becomes discontinuous, with the 
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FIG. 5: (Color) (a) Distribution functions p(A rn ) after integrations over long but finite times (t = 5 • 10 6 ) of initial conditions 
consisting of weakly perturbed 3D constant-amplitude (q x = q y = q z = 0) unstable solutions with a — 9 (blue), a = 12 (red) 
and a = 15 (black). System size 64 x 64 x 64 (N = 262144), a = C = 1. (b) Intensities, in a representative 15 x 15 x 15 subbox 
of the simulation box, at the end of the simulation for a — 9 in (a). Red and yellow patches are localized breathers. 



part below the breaking point corresponding to a phonon bath at T — oo, and the points above to large-amplitude 
breathers with increasing amplitudes. This is illustrated by Fig. ^d). Thus, this suggests that the separation of 
phase space into two parts as proposed in Ref. is valid also for larger a, although, as discussed in previous 
subsections, the time-scales to actually reach a true equilibrium state may be enormous and beyond reach of any 
numerical simulations. 

It should be obvious that also the extension to 3D is straightforward. We can e.g. consider a travelling plane wave 
in a cubic lattice, ?p mx ,m m, = y^ e '(fe"' I +«s" 1 !/+fe'»!) e !Af i an d obtain immediately the location of the localization 
transition line by adding the term cosg z to the numerator of l|16l) . Taking a — 1 and q x = q y = q z = 0, the critical 
value then becomes a = 12C for a constant-amplitude solution in 3D. This is illustrated numerically in Fig. [SJ Again 
we see that the distribution (Fig. EJa)) has the expected curvature both in the breather- forming regime (blue circles) 
where /3 < and in the normal regime (black circles) where (3 > 0. In Fig- Etb) we see that high amplitude breathers 
indeed do exist in the system for a = 9. 

To conclude this section, we thus see that, in contrast to the condition for existence of an energy threshold for 
creation of a single breather, which only involves the product aD, there is no equivalence between the spatial dimension 
and the degree of nonlincarity as concerns the existence of an equilibrium state with persistent breathers. Indeed, 
the presence or absence of such a threshold only affects the approach to equilibrium and not the qualitative features 
of the equilibrium state itself. The degree of nonlinearity and the dimensionality in our case actually tend to work 
in opposite directions, as we have seen e.g. for a constant-amplitude initial condition ip n = that increasing a 
decreases the maximum amplitude a for which persistent breathers form (see 1)12(1 ). while increasing the dimension 
increases it. 



III. KLEIN-GORDON CORRESPONDENCE TO DNLS PHASE TRANSITION LINE 



Let us now discuss how the DNLS statistical localization transition manifests itself for general KG chains of coupled 
classical anharmonic oscillators. In order to derive approximate expressions for quantities corresponding to the DNLS 
Hamiltonian and norm densities valid for small amplitudes and weak coupling, we follow the perturbative approach 
outlined in Ref. ^| (see also Ref. H2)- The KG Hamiltonian H for a chain of N oscillators is given by 



N r 



H 



E 

71=1 



-U 2 n + V(u n ) + -C K (u n+ l - U n ) 2 



where the general on-site potential V(u) for small-amplitude oscillations can be expanded as 

V(u) = -u 2 + a —+l3'— + .... 

The KG equations of motion then take the form 

tin + V'{u n ) - C K (u n+ i + U n _! - 2u n ) = 0. 



(17) 



(18) 



(19) 
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Considering small-amplitude solutions u n (t) with typical oscillation amplitudes \u r , 
panded in a Fourier series as 



in(t)=$>» )e 



e. they can be formally ex- 
(20) 



where ub is close to some linear oscillation frequency and the Fourier coefficients are slowly depending on time, 
a,n \e 2 t). Due to exponential decay of the Fourier coefficients in p they must satisfy ~ eP f° r P > 0, while 
an°' ) ~ e 2 - Moreover ai P ' = a„ since u„ is real. Inserting (|20|l into (|19fl yields: 



2 [a£> + 2ipco b aM + (l- P W b )a^ - C K {a^ +1 + a^ - 2<#>) 



a ipu b t 



+a 



(p) e ipu/ b t 







(p) e ipuj b t 



O + 0(e 4 ). 



(21) 



Then, we derive from I121H for the respective harmonics p = 0, 1, 2, the three equations ^ 



<4 0) + 


2«k«| 2 




-a (0) 
a n-l 


- 2«ff)) 


= o^ 


-C)(e 4 ), 


2z Wb d« + (1 - ^) a« + 2a(a^a^ + a^a^) + 3/3' 




- C K {a n \i - 


-a (1) 


- 2aW) 


= 0H 


-0(e 5 ), 


(1 - M 2 )a[f) - 


fa(aW) 2 


- C K (a% - 


r 0, n -l 


- 2«i 2 )) 


= 0- 


fO(£ 4 ) 



Consider first the case of a symmetric potential. Then, all odd powers of u in the expansion 118|) vanish (implying 
a = and C(e 5 ) in l|21|l ). and we immediately obtain a DNLS equation to 0(e 5 ) by considering Eq. I|23|) for the 
fundamental harmonic p = 1. 

For the general (non-symmetric) case, we proceed as in Ref. ^3j by assuming weak coupling Ck ~ e 2 (note that 
this assumption is not necessary to derive the DNLS equation for the symmetric case). Then, we can solve (I22|l to 
obtain: 



and l|24[) to obtain 



a^^-2a\a^\ 2 +0(e% 



ai 2 ) = |(aW) 2 + 0(e 4 ). 



(25) 



(26) 



(These are the weak-coupling limits of the more general solutions (15)-(18) in Ref. |l3|.) Inserting l|25|) - (|26|l into 
we get the general DNLS equation to 0(e 5 ) 

2iu b a^ + (1 - „ 2 )a« - C K (a% + - 2a«) + (-fa 2 + 3/3') l^^l 2 ^ 1 ) = + 0(e 5 ). (27) 

Defining 8' = c ^,~ 1 , A' = — ^a 2 + 3/3', er' = sign(X'), redefining time as t' = ^j-t, rescaling the amplitudes and 

moving into a rotating frame by defining ip' n — y^On'e 1 ^ ~ 2 ^* and neglecting terms 0(e 5 ), the DNLS equation in 
the new (slow) time variable t' takes the standard form 



(28) 



equivalent to with er = 1. For l|28|l we have the familiar conserved quantities as norm A — Y^, n =i IV4I 2 an d 



Hamiltonian H = J2n=i ( + "<Mn+\ ~ ^WS) ■ with h =< H > /N , a =< A > /N as before, the transition 
curve © between breather-forming and non-breathcr-forming regimes becomes h = —a' a 2 (breather-regime is above 
for a' = — 1 and below for a' = +1.) We now wish to express this condition in KG quantities. First, we express the 
norm as 
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By taking a£>* ■ $2$ - ■ l|?7)l* and summing over n, we find ^ \ T in |al 1} | 2 ) - e 6 iV, which together with Q 

2 

implies that the DNLS-norm in the general KG model behaves as A/N ~ ■^f(e t) (where / is some function of 
order 1). The DNLS Hamiltonian is then expressed as 



(30) 



By taking dl 1} * • (JUJl+d^ • O*, summing over n and defining ft« = E„[-2%n 3 | 2 + T I"" 5 1 4 + C ^ l a i+i ~ °»' I 2 ], 
where 5 = (w 2 — 1) /2, we find d7 ^| - ~ e 8 iV. Imposing the assumption of small coupling Cjc ~ e 2 (which together with 
the small- amplitude condition also implies 5 ~ e 2 ), we get that W/iV = — ^1 + 2(<5 — C^r) J2 n \ a n \ 2 ] ~ /(e 4 *)- 
Thus, the DNLS quantities .A/iV and W/iV correspond in the general case to two KG quantities of order unity, whose 
time variation is (at least!) two orders of magnitude slower than the typical time scale for the Fourier amplitudes 
(which in turn is two orders of magnitude slower than the time scale of oscillations of the original amplitudes u n ). 

Let us now explicitly calculate these quantities in terms of KG amplitudes and velocities u n , ii n . We do this by 
calculating time-averages of the different contributions to the KG Hamiltonian (|17(l with general potential energy 
(|18|l . Inserting the expansion H2U|) . averaging out all oscillating terms and using (|25|) - (|26[) . we get 

< E f >= \ E < f E a ^ bt + °^)) >= E (Vi 1 ^ 2 + \^ ] ) 2 + ^i 2 ) + o{j) 

n n \p=— 3 / n ^ 



= El^ 1) l 2 + y« 2 El a i 1) l 4 + ^ 6 )- (3i) 

n n 

Further, using also l |27f t we get for the time-averaged kinetic energy 

< E f >= \ E < f E <H P) + ^ bt + °(^) > 

n n \p=— 3 / 

= (\a^\ 2 + 4|4 2 ) | 2 ) + iu> b £ (<W - a«*d«) + 0( e 6 ) 

n n 

= (1 + 2C K ) J2 I4 1} | 2 + (-f " 2 + 3/3') ^ |a«| 4 - C K £(a« l0 £>' + + 0(e% (32) 

For the time- averaged cubic energy we get 

<£«y>=fE<(£# )e ^ t+0 ^>) > 

n n \p=— 3 / 

= f E (fttfWr' + 3 (4 2) («1 1} *) 2 + al 2) *(4 1} ) 2 )) + o(e 6 ) - - y« 2 E l^l 4 + °( £6 )' ( 33 ) 

n n 

for the quartic energy 

/ 3 \ 4 



<E/ 3 T >= tE< E + >= 2?' E Itf 5 1' + °( £6 )' (34) 

n n \p— — 3 / n 



and for the coupling-energy 



< E - «») 2 >= 2C '< E i^i 2 - ^ + + <^ 6 )- (35) 

n n n 

Using JHJ|, we can then write an approximate explicit expression for the DNLS norm as: 

^(<£f>+|<s4>W m 
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Note that in particular for the symmetric case (a = 0), the DNLS norm is, to 0(e 4 ), directly proportional to the 
averaged harmonic part of the on-site potential, while for the general case there is also an additional correction due 
to the cubic contribution. 

Then, there are several (indeed, infinitely many) ways of combining the quantities i|31|) - (|35[) . which all yield ap- 
proximate (to order e 2 ) expressions for the DNLS Hamiltonian i|30(). One way involving the KG Hamiltonian H i|17|) 
(showing that H and TL indeed are nontrivially related) is to write 



H 



H-<Y. ! T>-( 1 + 2C k 



Ewi 1 \ aui 



OU 2 



(37) 



This is in some sense the most appealing KG analog to the DNLS Hamiltonian, since it emphasizes the contributions 
from the coupling- and quartic energies to the KG Hamiltonian. Using this expression, we obtain the condition for 
the phase transition curve in terms of KG Hamiltonian and other quantities as: 



H 



II , I 1 T-s II. I-' I V ' II., \ I. Y '/„ 

— = A — < > — > H < > a— > H < > — 



19 1 



N 



N 



30 N 



N 



> 



+0 



(38) 



An example of another expression for TL is 
|A'| 



H = - 



2C% 



ff-2(l + 2C A -)<£^>-f <E 



au,. 
~3~ 



> 



<EY > 



(39) 



which notably does not explicitly include the quartic part of the on-site energy. 

Note also the following: By adding together all contributions from (|31() - (|35|l . we express the KG Hamiltonian H in 



terms of the fundamental Fourier amplitudes 

H = 2(1 + 2C K ) k X) | 2 + (-f o? + 9 -A Y k 1} | 4 - 2C K Y (oSi^ + + 0(e% 

n ^ ' n n 

Comparing with the expression l|3U|l for the DNLS Hamiltonian, we see that, generally, 



H 



3, 



2(1 + 2C K )J2\ 



,(i)| 2 
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/?'K:i a wi 4 + o( e 6 ). 



(40) 



(41) 



So in the very special case when a 2 = the coefficient in front of 



,(1)|4 



in 1)41(1 vanishes, and then (and only 



then!) is it possible to simply express the KG conserved quantity H in terms of the DNLS conserved quantities Ti 
and A: 



H = 



2C K 
|A'| 



(-C K H+(l + 2C K )A) + 0(e 6 ), 



(42) 



and to obtain an expression for the phase transition curve involving only the average KG Hamiltonian H/N and the 
average norm A/N (calculated e.g. using (JSSJ): 



H 

N 



(43) 

l) 2 , belongs to 



It is quite remarkable that one of the most studied examples, the Morse potential V(u) = \{e~ 
this special class, since for Morse a = —3/2 and j3' = 7/6 (=> A' = —4). 

In Fig. we show an example of results from long-time numerical integration of the Morse KG model, with a slightly 
perturbed constant-amplitude solution as initial condition. As is well-known, such an initial condition leads to breather 
formation through the modulational instability (e.g. Ref. H2), which is explicitly shown in Fig. EJa). In Fig. HJb) 
we show the variation of the above-derived approximate expressions for the DNLS quantities A and TL during the 
simulation time. Note that for moderate integration times (middle part of Fig.^b)) the three different expressions for 
Ji are close and agree well within the expected accuracy C(e 2 ). They also remain far from the localization transition 
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FIG. 6: (color online) Numerical integration of KG Morse chain with Ck = 0.005, TV = 200, and randomly perturbed constant- 
amplitude initial condition u n (0) = 0.05. (a) Time evolution of local energy density (note logarithmic time scale), (b) Main 
figure: Ti/N vs. A/N for the simulation in (a), with A calculated from l|36|l and, from top to bottom in the left part of the 
figure, Ti calculated from 1391 [red], 1371 [blue], and 1421 [green], respectively. Time runs from right to left (i.e. A/N decreases). 
Lowest curve is the localization transition line @. Inset in (b) shows the ratio of time-averaged cubic (1331 to quartic 1341 
energies versus time, compared to the DNLS prediction 1441 (lower line). 



line (lower curve in Fig. HJb)). However, for larger integration times the three curves diverge from each other 
(left part of Fig. EJb)), where in particular (|42[1 and <|37|) indicate an asymptotic decrease of Ti while l|39[) indicates 
an increase. This discrepancy can be traced to the fact, that the different expressions give different relative weights 
to the cubic and quartic anharmonic energies. As long as the amplitude remains small everywhere in the lattice, 
this difference is not important as all expressions are equivalent to C(e 2 ). However, as breathers grow, locally the 
oscillation amplitudes become significantly larger, indicating the beginning of a local breakdown of the validity of the 
DNLS approximation at the breather sites. According to H33|) - (I34[1 . the ratio between the averaged cubic and quartic 
parts of the anharmonic on-site energy remains fixed within the DNLS approximation, 

<e4 >/< ^ /3 t > =-?^ +o(£2) - (44) 

n n 

As can be seen from the inset in Fig. Elb) , the relative contribution from the quartic energy continuously increases 
with time, and gets significantly larger than the DNLS prediction (|44|l as the breathers grow. 

As another illustration of the role of the DNLS quantities for the KG dynamics, we consider a thermalized KG 

4 

lattice with a pure (hard) quartic potential, V(u) = ^- (i.e., a = 0, /3 = 1). We perform the following numerical 
experiment. First, we drive the system into a thermalized state by coupling it to a thermal bath at temperature 
T' , using standard Langevin dynamics by adding a fluctuation term —F n {t) and a damping term r\u n to the left- 
hand side of p9(l . (Note that this temperature T' is not equivalent to the previously discussed DNLS temperature 
T, since, as shown above, the DNLS Hamiltonian Ti is non-trivially related to the energy H of the KG-chain.) 
The fluctuation force F n (t) is taken as a Gaussian white noise with zero mean and the autocorrelation function 
< F n (t)F n > (t 1 ) >— 2r]T'd(t — t')S nn >, according to the fluctuation-dissipation theorem (with ks = 1). As can be 
seen from Fig. 0a), with the chosen damping constant r\ = 0.1 the lattice thermalizes after a few thousands of 

time units, with a time-averaged total kinetic energy < JJj ^ >= -jT 1 as expected. In the thermalized regime 
(t > 4000 in Fig. UJ, we monitor the quantities A/N and Ti/N calculated as instantaneous time-averages over fixed 
time-intervals (< f(t) >= ^ f t t f(t')dt' , where t = 100 in Fig. [7J. The results for a large number of time instants 
are illustrated by the dots in Fig. Efb). Note that taking simultaneously the limits (3'T' — > (harmonic oscillations) 
and Ck ~ * (thermalized uncoupled oscillators) with ^j^- constant, Eqs l|36() and (|37|l (or l|39|l ) yield ^ — > 

and ^ — > — | (%^-) ! which for the parameter values of Fig. |7Jb) corresponds to the point (0.75,-0.5625) on the 

localization transition line (dashed line in the figure). As can be seen, the effect of small but nonzero coupling and 
anharmonicity is to shift the long-time averages (center of the 'cloud' of dots in Fig. EJb)) towards smaller ^ and 
larger ^ (approximately (0.725,-0.505) in Fig. |7fb)), moving slightly into the 'non-breather-forming' regime of the 
DNLS approximation. However, due to the continuous interaction with the heat bath the fluctuations are large, and 
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time Norm density 

FIG. 7: Thermalization of a quartic KG chain (a = 0,f3' — 1) with C'k ~ 0.01, N = 800, coupled to a thermal bath at 

. 2 

temperature T' = 0.005 with dissipation constant r\ = 0.1. (a) Time-averaged total kinetic energy < ^p- >. (b) H/N vs. 
A/N for the simulation in (a), with A calculated from l|M6[l . and H calculated from IfiT). Each dot represents a time-average 
over the interval [t — 100, t] at 15382 different times t. Line in (b) is the localization transition line JSj. Larger points in (b) 
show the locations of the initial conditions used in Fig. |H1 
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FIG. 8: (a) Example of a breather appearing in the microcanonical integration of an initial condition represented by the lower 
large point in Fig. 0(b). (b),(c) (Non-normalized) velocity distribution functions p(u n ) obtained from long-time numerical 
microcanonical integrations (points) compared to Maxwellian distributions P(u) ~ (2?tT')~ 1 ' /2 exp(— u 2 /2T') (lines) at the 
estimated temperature. In (b) the initial condition is the same as for (a), the temperature is T' fa 0.00494 and the integration 
time is 1.2 • 10 6 . (c) corresponds to the upper large point at (0.704, —0.471) in Fig. 0b), with T 1 ~ 0.00485, and integration 
time 0.6 • 10 6 . In both cases, the velocities of all sites are registered in intervals of 0.6 time units. Insets in (b) and (c) show 
magnification of the small-velocity regime in non-logarithmic scale. 

the probability to be in the 'breather- forming' regime (below the dashed line in Fig. 0b)) at a given time-instant 
considerable. 

We then consider the effect of turning off the heat bath in the simulations in Fig. at different time instants, 
and continuing a microcanonical integration with the corresponding thermalized state as initial condition. We first 
choose an initial condition in the 'breather-forming' regime, corresponding to the point (0.716, —0.568) in Fig. 0b). 
(Even though this point is below the 'cloud' of dots in Fig. 0b) it does not represent a particularly exceptional initial 
condition in the thermal ensemble, since the dots represent time-averaged values rather than instantaneous, and the 

. 2 

fluctuations of the latter are considerably larger.) For this particular initial condition, monitoring < J2 n if - > during 
the microcanonical integration shows that it corresponds to a lattice temperature T' sw 0.00494. It is quite remarkable, 
that even with integration times longer than 10 6 we observe no systematic drift of either of the quantities A/N or 
H/N. Moreover, the fluctuations of these quantities calculated as fixed-interval time-averages over 100 time units as 
in Fig. 0(b) are very small (less than 5 • 10 -4 for A/N and 3 • 10 -3 for H/N ) and practically negligible on the scale 
of Fig. 0b). Thus, the system will remain in the 'breather- forming' regime, at least for extremely long time-scales. 
Although most of the breathers that can be observed are rather small and short-lived, examples of larger breathers 
persisting for about 20000 time units or more are not unusual and appear repeatedly throughout the integration time 
(see an example in Fig.|Sfa)). 
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To further illustrate the dynamics on the two sides of the transition line in Fig.[7|K), we compare in Fig. 0(b), (c) the 
velocity distribution functions p{u n ) obtained by long-time integration of two initial conditions corresponding to the 
two large points in Fig.Efb). In the breather-forming regime (Fig. Et'b)) the calculated p(u n ) shows a clear deviation 
from the standard Maxwell distribution, with a significantly enhanced probability of larger velocities (0.2 < \u n \ < 0.35 
in Fig. IHtb)). Also the probability of very small velocities {\u n \ ^ 0.04) is enhanced (see inset in Fig. 0b)), while the 
probability for intermediate velocities is decreased compared to the Maxwell distribution (the decrease for \u n \ > 0.35 
in Fig. 0b) is likely to be related to the finite size of the system) . Thus, the breather-forming processes tend to 
polarize the lattice into 'hotter' regions of larger oscillations and 'colder' regions of smaller oscillations, although due 
to the repeated creation and destruction of breathers at different sites, the equipartition result < >— T'/2 is 

still valid for each site, provided that the time-average is taken over a sufficiently large interval. On the other hand, 
for the initial condition belonging to the 'non-breather-forming' regime (Fig. 0c)), no such polarization relative to 
the Maxwell distribution can be observed. 



IV. CONCLUSIONS 

We have shown how a statistical-mechanics description of a general class of discrete nonlinear Schrodinger models 
yields explicit necessary conditions for formation of persistent localized modes, in terms of thermodynamic average 
values of the two conserved quantities TL and A. Furthermore, we illustrated how this approach can be extended to 
approximately describe situations with non-conserved but slowly varying quantities (see also Ref. [16J for a different 
example), and explicitly used it to explain formation of long-lived breathers from thermal equilibrium in weakly 
coupled Klein-Gordon oscillator chains. Concerning the roles of the degree of nonlinearity a and lattice dimension D 1 
we found that, in contrast to the condition for existence of an energy threshold for creation of a single breather, which 
involves only the product aD, a and D tend to work in opposite directions as concerns the statistical localization 
transition. The energy threshold affects only the approach to equilibrium and not the qualitative features of the 
equilibrium state. 

There are several directions in which we believe that this work should be continued. One important issue is 
to develop a quantitative theory determining the time-scales for approach to equilibrium in the breather-forming 
regime. As we have seen numerically, these time-scales may be extremely long, and naturally one may argue that 
the equilibrium states themselves are not physically relevant if they can only be reached after times of the order of 
t ~ 10 60 . Another important point regards, whether the hypothesis of separation of phase space in low-amplitude 
'fluctuations' and high-amplitude 'breathers' in the equilibrium state on the breather-forming side of the transition 
can be put on more rigorous grounds. Our numerical simulations are not completely conclusive in all the studied 
cases due to extremely long equilibration times, but give indications that this hypothesis could be valid also for large 
values of the norm density a. 

Finally, we stress the important connections to current experiments: Very recently, unambiguous experimental 
observations of discrete modulational instabilities have been reported, for an optical nonlinear array |36j |. as well 
as for a Bose-Einstein condensate in a moving optical lattice |37l ]. It will be very interesting to see, whether such 
experiments also can confirm the DNLS result that the final outcome of these instabilities depend, in a qualitative 
and quantitative manner, on the particular values of the Hamiltonian and norm densities (the latter represents power 
in the optical case and particle density in the Bose-Einstein context) as predicted here. 
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